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A recent experiment [Lanting et al, PRX, (2014)] claimed to provide evidence of up to 8-qubit 
entanglement in a D-Wave quantum annealing device. However, entanglement was measured using 
qubit tunneling spectroscopy, a technique that provides indirect access to the state of the system 
at intermediate times during the anneal by performing measurements at the end of the anneal 
with a probe qubit. In addition, an underlying assumption was that the quantum transverse-field 
Ising Hamiltonian, whose ground states are already highly entangled, is an appropriate model of the 
device, and not some other (possibly classical) model. This begs the question of whether alternative, 
classical or semiclassical models would be equally effective at predicting the observed spectrum and 
thermal state populations. To check this, we consider a recently proposed classical rotor model 
with classical Monte Carlo updates, which has been successfully employed in describing features 
of earlier experiments involving the device. We also consider simulated quantum annealing with 
quantum Monte Carlo updates, an algorithm that samples from the instantaneous Gibbs state of 
the device Hamiltonian. Finally, we use the quantum adiabatic master equation, which cannot 
be efficiently simulated classically, and which has previously been used to successfully capture the 
open system quantum dynamics of the device. We find that only the master equation is able to 
reproduce the features of the tunneling spectroscopy experiment, while both the classical rotor 
model and simulated quantum annealing fail to reproduce the experimental results. We argue that 
this bolsters the evidence for the reported entanglement. 


I. INTRODUCTION 

The D-Wave processors [1-3] are designed to be physi¬ 
cal quantum annealers [4] performing adiabatic evolution 
using programmable superconducting flux qubits. These 
devices have generated a substantial debate in the quan¬ 
tum computing community by laying claim to being the 
first large scale implementation of a quantum algorithm 
[5]. Much effort has been directed at answering the fun¬ 
damental question of whether the D-Wave processors ex¬ 
hibit sufficient “quantumness” to justify these claims. 

Independent verification of the quantumness of the D- 
Wave devices is challenging in part because of their black¬ 
box nature: the user interacts with the device by pre¬ 
senting it with an Ising model problem instance that is 
programmed as input, and receives a classical bit string 
representing the measured state of the qubits in the com¬ 
putational basis as output, at the end of the compu¬ 
tation (or quantum annealing run). This state is the 
device’s attempt at finding the ground state of the in¬ 
put Ising problem instance. This input-output interac¬ 
tion mode is clearly not amenable to the usual tests of 
quantunmess emphasizing non-locality [6]. Nevertheless, 
for sufficiently small problems (< 20 qubits), specific in¬ 
stances that emphasize quantum features of the evolution 
have been designed and found to show strong agreement 
only with quantum master equations, i.e., open quan- 
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turn system models, but not with classical models [7-9]. 
However, for a class of much larger (> 100 qubits) ran¬ 
dom Ising model problems, the device’s output exhibited 
strong correlations [10, 11] with a classical rotor model 
with Monte Carlo updates due to Shin, Smith, Srnolin 
& Vazirani (SSSV) [11], with simulated quantum an¬ 
nealing (SQA) implemented using quantum Monte Carlo 
[10, 12, 13], and with Parallel Tempering simulations [14], 
While a detailed study of the excited states and degener¬ 
ate ground states showed significant deviations between 
the device and these models [15], these results neverthe¬ 
less keep alive the question of whether SQA and the SSSV 
model provide an effective microscopic description of the 
device at large numbers of qubits. This question is par¬ 
ticularly pertinent since both models are efficiently sim- 
ulatable on classical computers. Another approach that 
has been used successfully to model the D-Wave devices 
is the quantum adiabatic master equation (ME) [16]. The 
latter is so far the only model that has successfully cap¬ 
tured all aspects of the “quantum signature” experiments 
reported in Refs. [7, 8]. Importantly, unlike SQA, this 
quantum model does not lend itself to an efficient clas¬ 
sical simulation. A related master equation (based on 
the noninteracting blip approximation [17]) was success¬ 
fully used to model collective tunneling in experiments 
involving a D-Wave device [9]. 

In contrast to the black-box approach, permitting ob¬ 
servations only at the end of each annealing run, recent 
experiments found evidence of entanglement generated 
during the evolution of the D-Wave devices from input 
to output, effectively opening the black box [18]. Specif- 
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ically, the experiments reported in Ref. [18] showed, us¬ 
ing qubit tunneling spectroscopy [19], that the measured 
quantum spectrum and thermal populations are in strong 
agreement with the quantum spectrum and Gibbs state 
of the transverse-field Ising Hamiltonian—the Hamilto¬ 
nian the device is supposed to evolve under. This, in 
turn, allowed Ref. [18] to demonstrate the existence of 
entanglement in the device using negativity [20, 21] and 
an entanglement witness [22-24]. 

However, an important caveat is that these experi¬ 
ments only provide an indirect way to detect entangle¬ 
ment, since the underlying assumption is that the quan¬ 
tum transverse-field Ising Hamiltonian is an appropri¬ 
ate model of the device. That is, entanglement was de¬ 
tected under the assumption that the measured spectrum 
and populations arise from transverse-field Ising models 
whose ground states are already highly entangled, and 
not from some other (possibly classical) model. This begs 
the question of whether alternative, classical or semiclas- 
sical models would be equally effective at predicting the 
observed spectrum and thermal state populations. For, 
if so, the deduced entanglement witness would not be 
applicable. 

To make this important point clearer, we might for¬ 
mulate it as a “model loophole,” in the tradition of the 
loopholes associated with nonlocality and Bell inequal¬ 
ity tests (see, e.g., Ref. [25]). The loophole is simply 
the fact that entanglement was detected under the as¬ 
sumption that a particular, already-quantum model of 
the dynamics, is responsible for the observed experimen¬ 
tal results. In an attempt to close this loophole, here 
we numerically simulate the tunneling spectroscopy ex¬ 
periment using the SSSV model, SQA, and the ME. We 
shall demonstrate that SQA and the SSSV model both 
fail to capture any feature of the entanglement witness 
experiments reported in Ref. [18], whereas the ME again 
succeeds. In other words, the reported measurements [18] 
are consistent with the ME, but inconsistent with SSSV 
and SQA. If the correct model of the D-Wave device is 
the ME then the reported measurements imply entan¬ 
glement. The SQA fails not because it isn’t a quantum 
model, but because it is the incorrect quantum model for 
the experiments at hand: unlike the ME it does not in¬ 
clude a unitary dynamics component, and this dynamics 
is what is presumably responsible for the experimental 
observations via an adiabatic connection of eigenstates. 
Our strategy does of course not rule out the possibil¬ 
ity that other classical models, or efficiently simulatable 
quantum models, might also obtain agreement with the 
reported measurements. Our results may be seen as an 
invitation to invent such models, which, if found, would 
further guide the study of the quantumness question of 
the D-Wave devices. 

The structure of this paper is as follows. In Section 
we review the principle behind the qubit tunneling spec¬ 
troscopy technique, as well as the theory behind the en¬ 
tanglement measures. In Section we describe the sim¬ 
ulation methods, in particular the ME, SQA, and the 


SSSV model. In Section we describe the simulated 
experiment. We present and discuss our results in Sec¬ 
tion V, where we demonstrate that only the ME matches 
the experimental tunneling spectroscopy results. We con¬ 
clude in Section . In the Appendix we demonstrate 
the robustness of our conclusions to various noise mod¬ 
els, and also reject another classical model based on spin 
dynamics with a friction term. 

II. REVIEW 

A. Qubit tunneling spectroscopy 

We briefly review the principle behind qubit tunneling 
spectroscopy [19], where the goal is to find the energy 
gaps of the quantum system Hamiltonian Hg. We take 
Hs to be of the form: 

N 

H s = -A^af + BH IS , (1) 

2=1 

where A, B > 0 are constants and 

His =Yl hia i + Y1 ( 2 ) 

i i<j 

is an Ising Hamiltonian acting only on the system qubits. 
The hi and Jtj are the local fields and couplings, respec¬ 
tively, and we use of (of) to denote the Pauli x (z) ma¬ 
trix acting on qubit i. 

We denote the eigenstates and eigenenergies of 
H s by {|-E„)}„=i and {£ n } ra=1 respectively with 
Hi < E 2 < • • •. A probe qubit P is coupled to 

system qubit 1 to give a system+probe Hamiltonian: 

H s+P =H S + BHip , (3a) 

Hip = Jipofop — Jipod — hpUp . (3b) 

An offset local field oc — Jip has been applied to qubit 1 
such that, in the eigenenergy subspace where the probe 
qubit is in the state |0), the eigenstates of the system are 
given by \E n ) ® |0) with energy E n — Bhp (where the first 
ket is the state of the system qubits and the second ket 
is the probe qubit state). When the probe is in the |1) 
state, the lowest energy state of the system and probe can 
be written as \ipo) ® |1) with eigenenergy eg = £o + Bhp , 
where |^o) is the ground state of Hs — 2BJppaf with 
eigenenergy eg- 

Let us assume that the system and probe are initial¬ 
ized in the state |'i/’o)® |1)- Introducing a small transverse 
field term (oc erg) for the probe qubit allows for transi¬ 
tions between the states \ipo) ® |1) and \E n ) (g> |0). In an 
open quantum system we may expect that the dominant 
process is incoherent tunneling between these two states 
[26]. By tuning the value of hp, we can make the two 
states degenerate, i.e., E n — Bhp = e 0 + Bhp, resulting 
in a resonant peak in the tunneling rate. Since both B 
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and hp are known, this allows us to solve for differences 
of the E n , and by finding the locations of the tunneling 
peaks as a function of hp , we can map out the quantum 
spectrum of Hg. For example, for a pair of such tunnel¬ 
ing peaks at Zip'* and hp\ corresponding to the n = 1 
and n = 2 energy eigenstates respectively, the energy gap 
between the eigenstates \E 2 ) and | E\) is then given by: 

E 2 -E 1 = 2 B (hp ] - 4 1} ) . (4) 


B. Equilibrium distribution 


Let us assume that to a very good approximation our 
system only populates the states |^o) <8> |1) and {| E n ) <g> 
|0)} n= i, such that the populations in these states sum to 
unity: 

Wo>®|l>) + £P(|E»>®|0>) = 1 ■ ( 5 ) 

n =1 


If we observe that the probe qubit is in the state |0), the 
system energy eigenstate populations P{E n ) are given 
by: 

p(F , = P(\E n )®\0)) = P(\E n ) <g> |0)) 

1 T,i=i p (\ E i) ® 1°)) 1--P(hfo>a|i» ■ 1 J 

However, if we have tuned hp so that the states \ipo) ® 11) 
and | E n ) (g> |0) are degenerate, and we have waited long 
enough that their populations have thermalized (i.e., 
they have equal populations), then P(\E n ) g> |0)) = 
P( IV’o) 8 |1)). Under these assumptions we can find the 
energy eigenstate populations entirely in terms of the 
population of the state \i/jo) <g) |1): 


P{E n ) = 


Wo)®|l)) 

l-P(\rPo)®\l)) • 


(7) 


We expect these to match the Gibbs state populations, 
i.e., P(E n ) = e~^ Hls /Z, where Z is the partition function 
and /3 = 1 /{kT) is the inverse temperature. 


C. Evidence for entanglement 

In Ref. [18], the authors use the populations found in 
the ground state Pi and first excited state P 2 to construct 
the density matrix of the system p = Pi\Ei)(Ei\, 

which assumes that the off-diagonal components in the 
energy eigenbasis are zero. Under this assumption 
(which, as we show below, agrees with the results of our 
ME simulations) the authors calculate the negativity [20] 
for all possible bipartitions A of the system, 

A f( P ) = \ (|| P Ta ||i -1) , (8) 

where p TA denotes the partial transpose of p with respect 
to the bipartition A. Global entanglement is then defined 


as the geometric mean of the negativity of all bipartitions 
[21], and was shown to be non-zero in Ref. [18]. 

A drawback of this approach is that it assumes the 
off-diagonal elements of the density matrix vanish. To 
ensure the robustness of an entanglement conclusion, an 
entanglement witness was used in Ref. [18] (based on the 
theory formulated in Ref. [24]): 

W A = \<t>){<t>\ TA , (9) 

where |</>) is the eigenstate of \E 0 )(E 0 \ rA with the most 
negative eigenvalue. The entanglement witness approach 
succeeds in certifying entanglement even when the off- 
diagonal elements of the density matrix are not con¬ 
strained to vanish, thus extending the range of validity 
of the presence of entanglement beyond that attainable 
using the negativity approach. 

Using this entanglement witness, two non-trivial 
checks were performed. First, experimental errors in the 
populations Pi and P 2 impose linear constraints on the 
state p: 


Pi - AP, < Tr [p\Ei)(Ei |] < P* + APi . (10) 

If Tr [Wap] < 0 for all p satisfying the experimental con¬ 
straints in Eq. (10), then entanglement is certified for the 
bipartition A. Optimizing Tr [Wap] subject to the above 
constraints is an instance of a semidefinite program, a 
class of convex optimization problems for which efficient 
algorithms are known. 

Second, uncertainties in the specification of the Hamil¬ 
tonian in Eq. (1), which can lead to changes in the eigen¬ 
states | Ei), were included by adding random perturba¬ 
tions (10 4 samples in total) to the Hamiltonian, while 
ensuring that the maximum of Tr [Wap] < 0 for all per¬ 
turbations. We independently verify that such noise does 
not change the presence-of-entanglement conclusion in 
Appendix 


III. SIMULATION METHODS 

The time-dependent Hamiltonian of the experiment is 
given by 


N 

H(s) = -As(s) ^2 erf - Ap(s)a x p + B(s)Hi sing , (11) 
2=1 

where s = t/tf is the dimensionless time, and the Ising 
Hamiltonian Hi slng is given by: 

Hising = His + Hip (12) 

Just as in Ref. [18], we take Jip = —1.8. The anneal¬ 
ing schedule functions (Ag(s), Ap(s), B(s)) are shown in 
Fig. . These schedules are calculated using rf-SQUID 
models with independently calibrated qubit parameters 
[27], and they correspond to the same schedule used in 
Ref. [18] [see their Fig. 1(d)]. 
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FIG. 1. (Color online) The annealing schedules used in the 
experiment [18]. The functional form for As(s) and Ap(s) in 
Eq. (11) is identical to the function A(s) shown. The dashed 
line corresponds to the experimental temperature of 12.5mK. 


A. ME 


strength gs to be the same and fixed for all the system 
qubits, and vary the probe-bath coupling strength gp rel¬ 
ative to gs , expecting gp > gs, since experimentally the 
probe qubit is operated in a regime where its coupling to 
the environment is strong (see Supplementary Informa¬ 
tion of Ref. [18]). In the simulations performed in this 
work, we fix the system-bath coupling strength to: 

gl'n/ri 2 = 1.2732 x 1CT 4 , (16) 

which is the value used in previous work [8]. 


B. SSSV 

The Hamiltonian of the SSSV model [11] is obtained by 
replacing erf i-a sin 6^ and erf i-a cosdi in Eq. ( .). The 
system is evolved by performing Monte Carlo updates 
on the angles £ [0,7r]. At the end of the evolution, 
the state is projected onto the computational basis by 
mapping Oi < 7r/2 —► +1 (state 0) and Oi > 7t/2 -a — 1 
(state 1). This model was derived from first principles 
using the Keldysh formalism in Ref. [28]. 


We assume that the qubit system is coupled to a bath 
of independent, thermal, Ohmic oscillators via a dephas¬ 
ing interaction. The adiabatic quantum master equation 
[16] can be used to describe the system evolution in the 
weak coupling limit, and the evolution of the density ma¬ 
trix is given by: 

d Z . r r / \ 1 

j t e = — h \n(t M 

N 2 

+ H ft [ L ab,a{t)p, K] + h -C- ,(13) 

cn—1 a,b 

where the index a runs over the qubits, the indices 
a,b = 1,..., 2 N run over the instantaneous energy eigen¬ 
values e a {t) of the system Hamiltonian H(t), with Bohr 
frequencies ujba = [£b(t) — e a {t)\/h. The Lindblad opera¬ 
tors are given by: 


Lab,a{t) = {e a (t)\(7a\£ b (t)){e a (t)\£b{t)) , (14) 

and the function T encodes the bath correlation function: 


r M =-j(u)+iS(u) 


(15a) 


irrje 


—M/u c 


1 — e~ 


duo' , 

& w 


UJ — UJ' 


(15b) 


where V is the principal value. The important free pa¬ 
rameters in the ME are the coupling strengths between 
the system and probe qubits to their respective bosonic 
baths. We denote these couplings by gs and gp respec¬ 
tively. For convenience, we take the system-bath coupling 


C. SQA 


We implement a discrete-time version of SQA [10, 12, 
13]. At each fixed time t in Eq. (11), Monte Carlo sam¬ 
pling is performed on the the dual classical spin system 
with: 

^ T ^ ] Jij 

i i<j 

^ , (1^) 

i.T 

where (3 is the inverse temperature of the Monte Carlo 
simulation, N T is the number of Trotter slices used along 
the Trotter direction, denotes the ith classical spin 
on the rth Trotter slice, and J T is the nearest-neighbor 
coupling strength of the 1th qubit along its Trotter direc¬ 
tion and is given by: 

ln(tanh (/3Ai(t)/N T )) > 0 . (18) 

In our simulations we fixed N T = 128 (we checked that 
increasing it does not alter the results). 


PHsV) = 


IV. DESCRIPTION OF THE SIMULATED 
EXPERIMENT 


The initial Hamiltonian is given by H( 1), with H(s) 
as in Eq. ( ). We choose the initial state of our ME and 

SQA simulations to be |1) = |1)®( JV+1 1 (the N system 
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hp h\ — — Jip h*2 — 0 



FIG. 2. (Color online) Depiction of the 2 + 1-qubit Ising 
problem. Qubits are displayed as (labeled) green disks, with 
their local field value above them. Couplings are shown as 
solid lines connecting qubits with their value above the line. 
Values are picked according to the experiment in Ref. [18], 
i.e., Js = —2.5, Jip = —1.8, and hp is varied. 


qubits plus the probe qubit); for SSSV we correspond¬ 
ingly choose all the initial angles as &i = 7r. We em¬ 
phasize that this is not necessarily the ground state of 
H( 1) = B(l)Hi slng . Note that by doing this we are able 
to skip the state preparation step performed in Ref. [18]. 
The simulation of the experiment then proceeds as fol¬ 
lows: 

1. B(s ) and Hg(s) are evolved backward from s = 1 
to s = s* in a time iq. The choice of s* determines 
the quantum Hamiltonian whose spectrum we wish 
to study. 

2. Hp(s) is evolved backward from s = 1 to s = sp in 
a time r\. 

3. The system evolves under the constant Hamil¬ 
tonian H = -H s (s*) J2ieS a i ~ Ap(s p )cr£ + 

B(s*)Hi s i ng for a “hold” time r. Note that this 
means that the values of A and B in Eq. (1) are 
given by H(s*) and B(s*) respectively. 

4. Hp is evolved forward from s = sp to s = 1 in a 
time T\. 

5. B(s) and Hg(s) are evolved forward from s = s* to 
s = 1 in a time t\. 

The state of the system qubits and probe qubit are then 
read. Since the states |1) and j^o) ® |1) are adiabatically 
connected energy eigenstates, measuring the population 
change in the |1) state indicates how much incoherent 
tunneling [26] (to the iso-energetic state \E n ) (g) |0)) has 
occurred at s *. By repeating the experiment for different 
values of the hold time r and recording the probability 
Pm ( T > hp, s*), of observing the |1) state at the end of the 
experiment, we can extract the tunneling rate T(/ip,s*) 
by fitting Pj^ (r, hp, s*) to the function a + be~ r( - hp ’ s *' >T . 
The experiments are repeated for a range of hp values in 
order to find the location of the peaks in T. 

The simulations were performed, as in the experiment 
of Ref. [18], on two system qubits plus one probe qubit, 
depicted in Fig. 2 , and eight system qubits plus one probe 
qubit, depicted in Fig. 3. 



FIG. 3. (Color online) Depiction of the 8 + 1-spin Ising prob¬ 
lem. The local fields are zero, and we take the couplings to 
be Js = —2.5. The probe qubit couples as in Fig. . 


V. RESULTS 

A. ME: 2 + 1 qubit results 

We first analyze the 2+1 qubit system example studied 
in Ref. [18], using the adiabatic quantum master equation 
described in Sec. . We perform the procedure out- 



FIG. 4. (Color online) Tunneling rate calculated using the 
ME, SSSV, and SQA for s* = 0.339 and grp = 10<?s compared 
against the experimental results in Ref. [18] [see the center fig¬ 
ure in their Fig. 8(a)], shifted to the left by 14.5GHz to align 
the first peak with the ME results (a constant term in the 
Hamiltonian present in the experimental setup but not in the 
ME accounts for this shift). For the experimental results, Ict 
error bars are shown. The two dominant peaks correspond to 
resonances with the ground state and the first excited state; 
the magnitude of the peak grows as <;p -+ gs (shown in Ap¬ 
pendix ), but the peak position remains the same regard¬ 
less of the value of gp. For SSSV, we show the results with 
T = 12.5 mK and n = 5. 10 6 runs were performed for each 
hp value, and the population of the |1) state was determined 
by counting the number of times it occurred in these runs. 
For SQA, we show the results with T = 12.5 mK and ti = 5. 
We performed 10® runs for each hp value, and the population 
of the |1) state was determined by counting the number of 
times it occurred in these runs. The tunneling rate for the 
ME is measured in gs -1 , while for SSSV and SQA it is in 
inverse sweeps. 
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(a) (b) 


FIG. 5. (Color online) ME results for the 2 + 1 qubit problem, (a) Energy levels calculated using the tunneling rate peaks via 
the ME with gp = 10gs- E on the vertical axis label is E 2 , E 3 , or £+ The solid curves correspond to the theoretical result 
from diagonalizing Hs- Error bars (only visible at a few points) represent 95% confidence interval for the fit of the mean for 
the Gaussian used to estimate the position of the peak in the tunneling rate (see Fig. ). (b) Population fraction as calculated 
using Eq. (7) via the ME with gp = 10gs, r = 1.5ms, and a temperature of 12.5mK. The solid curves correspond to calculating 
the Gibbs state associated with Hg. We do not show the population of the higher states since they are < 10 -5 . 


lined in Sec. with t± = 10/iS, and we give an example 
of the tunneling rate observed at a particular s* in Fig. 
and compare it to the experimental results of Ref. [18]. 
We estimate the position of the peak by fitting the data 
points around the peak with a Gaussian. We find excel¬ 
lent agreement with the energy spectrum as computed 
directly from diagonalizing the Hamiltonian, as shown in 
Fig 5(a). We then extract the energy eigenstate popula¬ 
tion distribution using Eq. ( ) at the values of hp corre¬ 
sponding to the peaks in the tunneling rate, and repro¬ 
duce the theoretical Gibbs state; see Fig. 5(b . Therefore, 
by reproducing the quantum spectrum and Gibbs pop¬ 
ulations, the ME is able to reproduce the spectroscopy 
results of the experiments performed in Ref. [18]. 

While the relative position of the dominant peaks cor¬ 
responding to the ground state and first excited state 
agree very well, we observe that the experimental re¬ 
sults are broadened considerably. It is shown in Ref. [18] 
that the experimental broadening is dominated by the 
linewidth of the probe qubit, which is claimed to be 
strongly affected by low frequency (1//) noise [29]. We 
attempted to reproduce the broadening using a variety 
of methods detailed in Appendix i, including the incor¬ 
poration of low frequency noise, but unfortunately were 
unable to do so, and this remains an open problem. Fur¬ 
thermore, the positions of the third and fourth peak cor¬ 
responding to the third and fourth excited states do not 
match the experimental result well, but this can be un¬ 
derstood as being due to a breakdown of the two-level 
approximation of the rf-SQUIDs (see the Supplementary 
Information of Ref. [18] and also Ref. [30]). 


B. SSSV: 2 + 1 qubit results 

We now follow the same procedure using the SSSV 
model. Because this model effectively operates in 
a strong system-bath coupling regime, it thermalizes 
rapidly, so we are forced to make T\ very small (we take 
each to be only 5 Monte Carlo sweeps, where a sweep 
is a complete update of all spins). Otherwise the SSSV 
model quickly forgets the initial state before it reaches 
s*, and also forgets the state it relaxed to at s* when it 
returns to s = 1. 

We show the tunneling rate calculated using the SSSV 
model in Fig. . Only a single peak is observed over the 
entire range of hp values studied, and this peak does not 
occur at the same position as any of the ME peaks. It 
is perhaps not surprising that the SSSV model does not 
reproduce multiple tunneling peaks, as its energy spec¬ 
trum is continuous. The single tunneling peak repre¬ 
sents the thermal configuration of the classical rotors at 
s* that maximally depopulates the |1) state. However, 
as the SQA results will demonstrate next, the failure is 
ultimately due to the absence of unitary dynamics that 
adiabatically connects energy eigenstates. 


C. SQA: 2+1 qubit results 

Unlike the SSSV model’s continuous energy spectrum, 
SQA’s spectrum is discrete. However, as we show in 
Fig. 1, SQA also fails to reproduce the experimental tun¬ 
neling signature on the 2 + 1 qubit problem and has a 
strong similarity to the SSSV model’s result (also shown 
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in Fig. ). We note that SQA and SSSV also correlated 
strongly on the 108 random Ising spin problem studied 
in Ref. [10], as further corroborated and explained in 
Ref. [15]. 

First, we check whether this failure is due to SQA 
somehow failing to capture the thermal expectation val¬ 
ues of observables. However, this is not the case: when 
we hold SQA at a constant Hamiltonian it correctly re¬ 
produced the thermal quantum Gibbs state populations 
for the computational states, as shown in Fig. 6. Instead, 
the present failure of SQA is rooted in the absence of uni¬ 
tary dynamics. Namely, there is no adiabatic connection 
between the initial |1) state and the state j^o) ® |1) at 
the point s*. In order for the arguments presented ear¬ 
lier to work, it is important that during the anneal from 
s = 1 to s* and sp and then back to s = 1, the evo¬ 
lution remain adiabatic with no loss of population from 
the energy eigenstates, in order for the population of the 
state jl) to accurately track the population of the state 
|^o) <8> |1)- For the ME, this is plausible since all other 
states with the probe qubit in the 11) state at s = 1 are at 
much larger energies, i.e., they correspond to much higher 
energy eigenstates. Therefore, in a simulation with a uni¬ 
tary dynamics component such as occurs in the ME, we 
do not expect (nor do we observe) any of these states 
to be populated such that the only state with the probe 
qubit in the |1) state is the |1) state. However, SQA 
lacks unitary dynamics, so at best it can provide an adi¬ 
abatic connection between the thermal state at s = 1 
and at s = s*, which cannot reproduce the experimen¬ 
tal tunneling peak signature. We also see that SQA does 
populate other states with the probe qubit in the 11) state 
besides the jl) state. (Note that if we choose to define 
the tunneling rate in terms of the populations of states 
with the probe qubit down, our SQA tunneling curves do 
not change.) 



ooo 001 010 011 100 101 no 111 

Computational state 

FIG. 6. (Color online) The population of the 8 computational 
states for SQA with N T = 128 Trotter slices at 1000 sweeps for 
10 6 repetitions, compared to their population in the quantum 
Gibbs state. Both are evaluated at s* = 0.339 and hp = 1.03, 
which is very close to the resonance condition for the states 
|-Ei) ® |0) and ji/'o) ® |1). 


D. ME: 8+1 qubit results 

We extend our analysis to the 8 + 1 qubit problem de¬ 
picted in Fig. 3, which was also studied in Ref. [18]. We 
follow the same method as in the 2 + 1 qubit case, the 
only difference in our numerical simulations being that 
we truncate the energy spectrum to the lowest 16 energy 
eigenstates in order to reduce the computational effort. 
The results are shown in Fig. . We continue to find 
no agreement between the ME and SQA or SSSV. The 
experimental results exhibit two broad tunneling peaks 
whose centers match the ME results after a constant shift 
(as in the 2 + 1 qubits case). While the amount of broad¬ 
ening is the same in the 8+1 and 2 + 1 qubit experimental 
results, the limited range of hp values used (extending 
only up to 15GHz as opposed to 45GHz in Fig. 4) makes 
the agreement appear less impressive than in the 2 + 1 
qubits case. 



FIG. 7. (Color online) Tunneling rate calculated using the 
ME, SSSV, and SQA for s* = 0.284 as well as the DW ex¬ 
perimental results (shifted to align with the first ME peak) 
on the 8+1 qubit problem. For the experimental results, Ict 
error bars are shown. We show the results with T = 12.5 mK 
for all simulation methods. Otherwise, the parameters for all 
simulation methods are as for the 2 + 1 qubit case. The tun¬ 
neling rate for the ME is measured in /rs -1 , while for SSSV 
and SQA it is in inverse sweeps. 

The energy spectrum and the populations as derived 
from the ME results are shown in Fig. . The agreement 
between the ME and the results from diagonalization [in 
panel (a)] and the Gibbs state [in panel (b)] is again 
excellent. 


VI. CONCLUSIONS 

In this work we set out to check whether the conclu¬ 
sion that entanglement is present during the evolution 
of 2- and 8-qubit experiments using the D-Wave device 
can be explained using a classical rotor model (SSSV), 




























FIG. 8. (Color online) Gap and populations ME vs Gibbs state results for the 8+1 qubit problem, (a) The gap to the 
first excited state calculated using the tunneling rate peaks via the ME with gp = 10gs- The solid curves correspond to the 
theoretical result from diagonalizing Hs- Error bars (barely visible) represent 95% confidence interval for the fit of the mean for 
the Gaussian used to estimate the position of the peak in the tunneling rate (see Fig. .) (b) Population fraction as calculated 
using Eq. (7) via the ME with gp = 10<?s, r = 1.5ms, and T = 12.5mK. The solid curves correspond to calculating the Gibbs 
state associated with Hs- We do not show the population of the higher states since they are < 10 -5 . In our simulations, we 
truncated the energy spectrum to the lowest 16 energy eigenstates to keep the computational time tractable. We kept track 
of the populations along the evolution in order to make sure that this truncation does not lead to a loss of population. We 
checked that the results do not change when additional energy levels are included. 


SQA, or the quantum adiabatic master equation (ME). 
This question is pertinent since the evidence for entan¬ 
glement is indirect and assumes a quantum Hamiltonian. 
We found that neither the SSSV model nor SQA can 
explain the evidence. The failure of both models can ul¬ 
timately be attributed to the fact that neither accurately 
captures the adiabatic connection between energy eigen¬ 
states during the annealing evolution. In contrast, the 
ME reproduces the experimental tunneling spectroscopy 
signature, suggesting that the underlying transverse field 
Ising mode Hamiltonian, with its quantized energy spec¬ 
trum, is an appropriate description, along with the weak 
coupling approximation used to derive the master equa¬ 
tion, at least for the system sizes of up to 8 qubits we 
have considered. 

Our results support the presence of entanglement in 
the evolution of the D-Wave device, as concluded in 
Ref. [18]. Furthermore, the failure of SQA and the SSSV 
model to capture this result indicates the importance of 
real-time open system quantum dynamics in supporting 
this conclusion. We emphasize that it is important to 
have both the right thermalization as well as the right 
dynamics to reproduce the experimental results; for ex¬ 
ample, we have checked that a Langevin-0(3) model as 
used in Ref. [8] fails to reproduce the experimental results 
as well (see Appendix ). 

Nevertheless, it is important to note that while it 
becomes computationally prohibitive to use the ME to 
study systems larger than about 15 qubits, SSSV and 
SQA do not have this drawback and seem to capture cer¬ 


tain experimental results (in particular the ground state 
population) at scales of > 100 qubits [10, 11]. This ten¬ 
sion between the small and large system behaviors re¬ 
mains an open question: does the open quantum sys¬ 
tem description of the device remain valid at larger sizes, 
and is the agreement of SQA and the SSSV model with 
ground state populations due to the fact that they and 
the open quantum system description yield similar final- 
time statistics? If so, this could be an artifact of the 
problem instances studied so far [31], or the dominance 
of thermal noise in the annealing evolution. Or does the 
weak-coupling approximation actually break down, so 
that semi-classical descriptions such as the SSSV model 
become valid for the dynamics as well? Further progress 
in open quantum system modeling and more appropriate 
benchmarking tests are required in order to address these 
questions. 
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Appendix A: Robustness of the entanglement 

We consider adding Gaussian noise to all the fields and 
couplings in the 2-qubit system Hamiltonian [Eq. ('.)], 
i.e., we include an additional Hamiltonian of the form: 

-ffNoise(s) = A(s) (aiCTi + a 2 al) 

+B(s) (a 3 of + "b oi^cr^ er^), (Al) 

where a t ~ 7V(0,0.1). This choice is well above the typ¬ 
ical a = 0.05 ICE associated with the D-Wave Two pro¬ 
cessors. We show in Fig. 9 that for the 10 5 noise samples 
generated, the Gibbs state never has zero negativity [as 
defined in Eq. ( )]. Therefore, such a noise model is 
unlikely to invalidate the presence-of-entanglement con¬ 
clusion of Ref. [18]. 



Negativity 


FIG. 9. (Color online) Histogram of the negativity of the 
Gibbs state associated with the 2-qubit Hamiltonian of Eq. (1) 
when a noise term of the form given in Eq. (Al) is included. 
The red vertical line at 0.211 is the negativity for the noiseless 
case. A total of 10 5 noise instances are shown. 


Appendix B: Extended noise models 

While the ME is able to reproduce the positions of the 
tunneling rate peaks (see Figs, and 7), the experimental 
results show a broadening of the peaks that is absent in 
the ME results, indicating that an important noise source 
is missing from the ME simulations. In the Supplemen¬ 
tary Information of Ref. [18], it is made clear that the 
line shape of the first peak (corresponding to the ground 
state) in the two-qubit and eight-qubit systems are al¬ 
most identical to the probe qubit’s macroscopic resonant 
tunneling (MRT) profile. This suggests that the broad¬ 
ening observed is from noise on the probe qubit and not 
the multi-qubit system. In our simulations, we have at¬ 
tempted to model this by increasing the relative system- 
bath coupling strength between the probe and system 
qubits. We show in Fig. 10 that the broadening induced 
by this method is not sufficient to explain the observed 
broadening. 



FIG. 10. (Color online) Tunneling rate for the 2 + 1 qubit 
problem calculated using the ME for s* — 0.339 and variable 
gp/gs compared against the experimental results in Ref. [18]. 



2B(s*)h P /h (GHz) 

FIG. 11. (Color online) Tunneling rate for the 2 + 1 qubit 
problem calculated using the ME for s* = 0.339 and gp = 
10gs with a = 0.05 and a = 0.1 ICE on the probe qubit only 
and without ICE compared against the experimental results 
in Ref. [18]. 

We attempt therefore to reproduce this broadening us¬ 
ing several other possible noise models. Our noise model 
choices are phenomenological, and are designed to test 
whether such modifications are sufficient to broaden the 
tunneling peaks. First, we consider introducing Gaussian 
noise on the Ising local fields and couplings: 

Jij —> Jij + Af(0, cr) , hi —> hi + A/"(0, a ) (Bl) 

This is a relevant source of noise for the D-Wave proces¬ 
sors, often referred to as internal control error (ICE) [32], 
We run the ME with 1000 noise realizations for each ap¬ 
plied (ideal) bias h p . We then average the observed pop¬ 
ulation in the |1) state over these 1000 realizations and 
then extract the tunneling rate associated with the ap¬ 
plied (ideal) bias h p . We first consider the case where the 
noise is purely on the local field and coupling of the probe 
qubit. We show in Fig. results with a = 0.05 and 
a = 0.1, well above the typical a = 0.05 ICE associated 
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FIG. 12. (Color online) Tunneling rate for the 2 + 1 qubit 
problem calculated using the ME for s* = 0.339 and grp = 
10gs with a = 0.05 and a = 0.1 ICE on all qubits and without 
ICE compared against the experimental results in Ref. [18]. 



FIG. 13. (Color online) Tunneling rate for the 2 + 1-qubit 
problem calculated using the ME for s* = 0.339 and gp = 
10gs averaged over a [see Eq. (1 )] compared against the 
experimental results in Ref. [18]. 



FIG. 14. (Color online) The modifications due to the addition 
of telegraph noise [Eq. (E )] to the purely Ohmic spectral 
density function 7 ( 0 ;). 



2B(s*)h P /h (GHz) 

FIG. 15. (Color online) Tunneling rate for the 3 qubit prob¬ 
lem calculated using the ME for s* = 0.339 and grp = 10gs 
using purely Ohmic and Ohmic plus telegraph spectral noise 
compared against the experimental results in Ref. [18]. 


with the D-Wave Two “Vesuvius” processors. While the 
introduction of ICE in the simulations reduces the height 
of the tunneling-rate peaks and broadens them slightly, 
it is insufficient to account for the broadening observed 
in the experiment. Similar results are achieved if ICE 
is introduced on all qubits (see Fig. 2), suggesting this 
noise model is not sufficient to explain the broadening 
obesrved. 

A more restrictive (and less physically reasonable) 
noise model is: 


hi —y hi + ol , % — 1,2, (B2) 

where we average over three a £ {— 0 . 1 , 0 , 0 . 1 } values. 
The result is shown in Fig. . While this does not suffi¬ 
ciently broaden the peaks, it does cause the appearance 
of peak splitting, which is a feature of the experimental 
results. 

Finally, we modify the spectral density +(w) in Eq. (1 ) 


to include a low frequency component. We consider a 
form of telegraph noise [33], which we model as 


7tel(w) 


U! 1 

1 — e~P u uj 2 + 


(B3) 


in order to satisfy the KMS condition [34]. We set 
w m = 0.01,0.05 GHz 2 (this choice is motivated by nu¬ 
merical stability; making wir and wir smaller results in a 
significant slowdown of our simulations), and assume that 
this contribution to 7 (oj) has the same coupling strength 
g as the Ohmic component. We show in Fig. how this 
modifies the purely Ohmic case. As we show in Fig. 15, 
these modifications are counterproductive and act to nar¬ 
row the peak. We also included 1/f noise [with a spec¬ 
tral density of the form 7l// ( w ) = but 

this did not reproduce the experimental broadening ei¬ 
ther (not shown). 
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(a) (b) 

FIG. 16. (Color online) Results for the spin Langevin model described in Eq. ( 1) for the 3 qubit problem. Simulation 
parameters are \ = 10~ 3 and T = 12.5 mK. 1000 runs were performed at each h-p value. (Mp) denotes the average over these 
1000 runs for the 2 -component of the magnetization of the penalty qubit. 


Since we were unable to reproduce the observed broad- and 
ening we must conclude that our noise model is lacking 
a relevant feature, likely related to the breakdown of the 
weak coupling limit. Alternative methods (such as the 
non-interacting blip approximation [9], though it is re¬ 
stricted to two-level systems) may thus need to be em¬ 
ployed or developed to capture the broadening aspect of 
the experiment. 


Hi = 2 Ai(t)x - 2 B(t) ( hi + J2 J ‘ 


ZjMj 


(C3) 


Appendix C: Spin-Langevin Model 

Given the importance of the unitary dynamics compo¬ 
nent in reproducing the experimental results, we consider 
an alternative classical model with dynamics, specifically 
a (Markovian) spin-Langevin equation [35, 36] with a 
Landau-Lifshitz friction term [36, 37], 

= ~(H t + at) + X Hi X Mi) x Mi , (Cl) 

with the Gaussian noise £ = {£i} satisfying 
<&(*)} = 0 , <6(t)&(i , )> = 2 k B T x 5 l3 5(t - t') , (C2) 


where x and z are unit vectors. Like the SSSV model, 
this model is another classical limit that can be derived 
from the Keldysh path integral formalism [28] . We follow 
the same procedure as described in Sec. of the main 
text. We find that magnetization along the 2 -direction 
of the penalty qubit does not depend on the hold time 
r, as shown in Fig. . Therefore, we do not observe 
the experimental signature of exponential decay of the 
all-down state. Furthermore, the final outcome of the z- 
magnetization in fact depends smoothly on hp as shown 
in Fig. . Therefore, although this model does in¬ 
clude dynamics, it fails to reproduce the experimental 
signature. 









